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I. INTRODUCTION 

The description of strongly correlated bosonic systems 
is of fundamental interest in largely diverse physical sit- 
uations ranging from low-temperature experiments with 
superfluid helium |l| to Josephson-junction arrays [|| , as 
well as magnetic insulators |3j and ultracold gases in op- 
tical lattices H, Q • The latter systems offer an unpar- 
alleled playground to study fundamental models widely 
considered in statistical and condensed-matter physics. 
This is because of the high degree of control over the ex- 
perimental parameters that determine the Hamiltonian 
describingthe system. In particular, the Bose-Hubbard 
model [flj3 bas been experimentally realized in one Q, 
two 0, [H|, and three dimensions where the su- 

perfluid to Mott-insulator transition has been observed. 
Even though it has received less attention, the super- 
fluid to normal transition in the Bose-Hubbard model 
has been investigated experimentally in three dimensions 
fl2j . while in two dimensions it has been realized in the 
form of a two-dimensional lattice of Josephson-couplcd 
Bose-Einstein condensates [TH, [l4| , as well as in experi- 
ments with ultracold atoms in optical lattices (l5j . 

Although experiments with ultracold atoms on opti- 
cal lattices are in some respects almost ideal realizations 
of model Hamiltonians of interest, significant complica- 
tions arise because of the presence of a confining poten- 
tial, which leads to the coexistence of different phases 
in a single experimental setup fl6L f]~7j . Furthermore, the 
mesoscopic size of the system in combination with the in- 
homogencity induced by the trapping potential produces 
a rounding off of the otherwise sharp features present in 
an infinite homogeneous system in the critical region 
|2~0 | . Thus the understanding and assessment of criticality 
in such systems remains a challenging task. 

The emergence of sharp features in the momentum dis- 
tribution as obtained from time-of-flight images has been 
frequently associated to the emergence of superfluidity 
[nHH 25]. However, this association may not be accu- 
rate because sharp peaks in the momentum distribution 



already appear in the normal state, due to an increas- 
ing correlation length when approaching a critical regime 
[2a - t28| . More recently, new schemes to detect criticality 
in trapped systems have been proposed. In some of those 
studies, a detailed analysis of the momentum distribution 
was used to define criteria that allow one to extract reli- 
able estimations of the critical points from time-of-flight 
images H3, ■ hi addition to time-of-flight images, 
high-resolution in situ imaging of the density profile of 
trapped systems has become a powerful instrument with 
which one can also study phase diagrams of strongly cor- 
related systems and quantum criticality. Numerous the- 
oretical and experimental studies based on this idea have 
been carried out for systems in the presence of an optical 
lattice [33443 and in the absence of it I41H43. 



One important aspect that determines the nature of 
the quantum phases and their associated order param- 
eters is the dimensionality d. Mermin et al. rigorously 
proved that at any nonzero temperature, continuous sym- 
metries cannot be spontaneously broken in systems with 
sufficiently short-range interactions in dimensions d < 2 
H^,|4^|. This implies that, at finite temperature, Bose- 
Einstein condensation (BEC) cannot occur in one and 
two dimensions. Two-dimensional Bosc systems, how- 
ever, are marginal in the sense that fluctuations are 
strong enough to destroy the fully ordered state but are 
not so strong as to suppress superfluidity. Thus critical 
behavior develops in the Berezinskii-Kosterlitz-Thouless 
(BKT) transition [13, |48| , where a superfluid phase with 
quasi-long-range order competes with thermal fluctua- 
tions and induces a continuous phase transition to the 
normal fluid as the temperature is increased. In addi- 
tion to low-temperature superfluidity, long-range order 
can develop at zero temperature in two dimensions. On 
the other hand, in three dimensions, the superfluid tran- 
sition is accompanied by the appearance of true long- 
range order, implying that the system also exhibits Bose- 
Einstein condensation. Such a transition, which belongs 
to the three-dimensional XY universality class, is well 
understood in the sense that the critical exponents have 
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been determined experimentally and theoretically with 
remarkably high accuracy in many different physical con- 
texts [HUE!- 

Here, we focus our study on the supcrfluid to normal 
transition in a system of strongly interacting bosons in 
two- and three-dimensional lattices. Specifically, we con- 
sider the Bose-Hubbard model in the limit of infinite on- 
site repulsion, i.e., the hard-core boson limit. We use 
exact quantum Monte Carlo simulations to compute the 
finite-temperature phase diagram as a function of chem- 
ical potential. Accurate results are obtained through 
finite-size scaling of the condensate fraction and/or the 
supcrfluid stiffness obtained from our simulations. We 
also determine the mean-field phase diagram, which is 
qualitatively correct but quantitatively quite different 
from the exact results. We then proceed to study the 
supcrfluid to normal phase transition in two and three 
dimensions in the presence of a confining potential, which 
is required to describe experiments with ultracold gases. 
We introduce a method to determine the critical temper- 
ature, for any given density, that is based on the mea- 
surement of the zero-momentum occupation as a function 
of temperature. This method is in principle adequate for 
experiments dealing with both homogeneous and trapped 
systems. Furthermore, we compare our approach to other 
recently proposed schemes based on the in situ density 
images [3l| as well as on the shape of the low-momentum 
part of the momentum distribution f29j . 

The paper is organized as follows. In Sec. UH we in- 
troduce the model and its phase diagram in two and 
three dimensions supplemented with the mean-field cal- 
culations. Section Hill is devoted to the discussion of the 
techniques to obtain the phase boundaries. In Sec. IIV| 
we discuss the possibility to have Bose-Einstein conden- 
sation in trapped two-dimensional systems as well as the 
methods to determine the phase boundaries from exper- 
imentally accessible quantities. Finally, in Sec. [Vj we 
draw our conclusions. 



II. MODEL AND PHASE DIAGRAM 

We consider a system of hard-core bosons on a d- 
dimensional lattice with L d sites. The Hamiltonian can 
be written as 

H = -t ^ (otoj- + H.c.) - ^ mfk , (1) 

where a\ (a, ) is the boson creation (annihilation) opera- 
tor at a given site i, and hi = a\a i is the local particle 
number operator. The hard-core boson creation and an- 
nihilation operators satisfy the constraint a] 2 = af = 0, 
which forbids multiple occupancy of lattice sites. The 
first term in Eq. ([T]) is the kinetic energy, where t is 
the hopping amplitude between neighboring sites i and j 
In experiments involving ultracold gases, a trap 
is required to confine the atoms. The effect is taken into 
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FIG. 1. (Color online) Finite temperature phase diagram in 
two and three dimensions, and the mean-field (MF) predic- 
tion. In all dimensions, the phase diagram contains a super- 
fluid (SF) lobe surrounded by the normal fluid (NF) phase. 

account in the second term that contains /i; = fi — VqV 2 , 
where Vo is its strength and (i is the overall chemical po- 
tential, n is the distance from site i to the center of the 
trap. In what follows, positions will be given in units of 
the lattice spacing a and the energy will be given in units 
of the hopping amplitude t. 

We recall that the Hamiltonian in Eq. ((I) can be ex- 
actly mapped to the extensively studied quantum XY 
model [561 ] 

H = -2tJ2 (sfsj + sfsy) ^>,.s; . (2) 

where Sf is the ath component of the spin-1/2 spin oper- 
ator at site i. In the spin language, the term proportional 
to t describes a ferromagnetic exchange interaction, while 
the one proportional to /i, describes a magnetic field in 
the z-direction at site i. 

We study the Hamiltonian in Eq. (p}, at finite tem- 
perature T, by means of the stochastic series expan- 
sion (SSE) quantum Monte Carlo (QMC) method with 
operator- loop updates (5714591 ]. The determination of 
the phase diagrams is carried out through a finite size 
scaling of the condensate fraction and/or the super- 
fluid stiffness p s using periodic boundary conditions. 
The numerically exact (QMC) phase diagram in two 
dimensions (2D) and three dimensions, as well as the 
the mean- field predictions, are presented in Fig. [T] 
The finite-temperature phase diagram comprises an off- 
diagonal long-range-ordered (ODLRO) low-temperature 
superfluid lobe (quasi-ODLRO in 2D) surrounded by a 
high-temperature normal phase with exponentially de- 
caying correlation functions. The extend of the super- 
fluid state is expected to be hindered as dimensionality is 
reduced because thermal and quantum fluctuations have 
a stronger effect in low-dimensional systems. Clearly, 
our results agree with that expectation. The dissimilar- 
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ity between the mean-field and the exact phase diagrams 
makes it clear that both thermal and quantum fluctua- 
tions are strong and play an important role even in three 
dimensions, where mean-field approaches are generally 
considered to be a good approximation. 

Details on the procedure to obtain the phase bound- 
aries are provided in the following sections. Such proce- 
dures are different in two and three dimensions because 
of the different universality class of the phase transition. 



III. HOMOGENEOUS SYSTEMS 
A. Two dimensions 

Our results for the two-dimensional phase diagram in 
Fig. Q] are based on the fact that the model in Eq. ([T]) 
undergoes a BKT transition as a function of the temper- 
ature. This phase transition has been studied in great 
detail the context of the two-dimensional quantum XY 
model in Eq. ((2J in the absence of a magnetic field (60l — 
l63l |. Kosterlitz and Thouless predicted that the super- 
fluid stiffness p s jumps from zero to the value (2/tt)T c 
at the critical temperature. Thus we consider measure- 
ments of the superfluid stiffness p s for different system 
sizes L as a function of temperature. Within the SSE 
method, the superfluid stiffness is computed by measur- 
ing the fluctuation of the winding number W [64| ; they 
are connected through the relation p s = (W 2 ) /2/3, where 
P = 1 /T is the inverse temperature. 

Figure [2ja) shows results for the superfluid stiffness of 
2D hard-core bosons at p = [or, equivalently, the spin 
stiffness of the 2D XY model in Eq. ([2])] as a function of 
T for several system sizes. The observed slow approach 
of the superfluid stiffness to the characteristic jump ex- 
pected for the infinite system is due to strong finite-size 
effects at the BKT transition. Finite-size scaling rela- 
tions for the superfluid stiffness can be derived by in- 
tegrating the Kosterlitz renormalization-group equations 
(see, for instance, Rcfs. 63.65.66). This procedure yields 
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p s (T c ,L)tt _ 

2T 2 (In L + l 

p s (T,L)ir 



1 = ccoth2c(lnL + Z ), T < T c 
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where c measures the distance from the critical point and 
lo depends only weakly on temperature. Close to the crit- 
ical point, c ~ y/\T-T c \. In the limit 2c (ln£ + 1 ) < 1, 
a scaling form for the superfluid stiffness based on Eq. ((3]) 
can be written as 
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From Eq. © in the limit 2c(lnL + Z ) < 1, F(x) = 
1 — (4/3) x. From Eq. one can find the scaling 
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FIG. 2. (Color online) (a) Superfluid stiffness in 2D for p — 
and several values of L. The error bars (not shown) are 
smaller than the point size used in the plot, (b) Data collapse 
according to the relation in Eq. @. The inset in (b) shows 
the rescaled superfluid stiffness vs T. 



function F and critical temperature T c by computing 
x L = (In L + lo f (T - T c ) jt and y L = p s (T, L) n/T - 2 
based on our Monte Carlo simulations for different L and 
T. The adjustment of the constant Iq and critical tem- 
perature T c , such that the data produce the best possi- 
ble collapse, yields a numerical estimate of the scaling 
function F and the critical temperature itself. The re- 
sult of the determination of the scaling function F is 
reported in Fig. [U[b), where a plot of ul as a function 
of xl is presented. Notice that as expected, the value 
of F is very close to one for xl = 0. Furthermore, one 
expects from Eq. ^ that a plot of the rescaled super- 
fluid stiffness p s (T, L)* = p s (T, L) (l + 2[ln £ +io] ) as 
a function of the temperature T should become system- 
size independent at the critical temperature T c . This 
observation is confirmed in the inset of Fig. HJb). Re- 
markably, those curves intersect with the line (2/ir)T 
right at the critical temperature, in agreement with the 
BKT scenario. Our result T c /t = 0.685 ± 0.001 is con- 
sistent with the best value reported in Ref. [6^, for which 
T c /t = 0.6846 ± 0.0006 [fi^. An analogous procedure 
to the one just described is carried out for different 
values of the chemical potential to complete the two- 
dimensional phase diagram in Fig.[TJ Wc should mention 
that Eq. predicts the value of the superfluid stiffness 
in an infinite system at the critical temperature to be 
p s (T c ) /T c = 2/tt. However, in Ref. Ir38l . it was shown 
that the superfluid stiffness at the transition tempera- 
ture is p s (T c ) /T c ~ 0.63650, which is very close to the 
result based on Eq. © (2/tt ~ 0.63662). Detecting the 
difference is beyond the accuracy of the present study. 



1. Critical value from dno/dT 

We now briefly discuss the behavior of the occupation 
of the zero momentum state (n^o = no) in the crit- 
ical region and address the determination of the tran- 
sition temperature from it. In a homogeneous and in- 
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finite 3D system, BEC is identified by a macroscopic 
occupation of n . However, as mentioned before, ther- 
mal fluctuations in 2D destroy Bose-Einstein condensa- 
tion. Nonetheless, as the superfluid transition is ap- 
proached from the normal phase, no diverges [see in- 
set in Fig. Eta)]. Indeed, from the Fourier transform of 
the one-body density matrix in the long-distance limit 
(a\a i+r ) oc r -1 / 4 exp (— r/£), one can extract the behav- 
ior of no as T c is approached, 



n 



£ 7/4 . 



(5) 



We assume the essential singularity of the correlation 
length £ ~ e b /VT-T c ^ waere 5 i s a chemical-potential- 
dependent scaling factor. From Eq. ((5J), it follows that 
not only does no diverge at T c , but also its derivative 
with respect to T, 



dno 
~dT 



c ? / 4 in 3 e 

b 2 



(6) 



In a finite system, when T is close to T c , the role of the 
correlation length is taken over by L when £ > L. This 
occurs at a characteristic temperature T* (L) given by 



T* (L) = T C + h' I In 2 L, 



(7) 



where b' is a nonuniversal factor related to b. At that 
temperature, the derivative in Eq. ^ scales with the 
system size as 



dno 
~dT 



T*{L) 



L 7 / 4 ln 3 L 
b 2 



(8) 



Below T*(L), no cannot vary as fast as right above 
T*(L) because the exponential increase of the correla- 
tion length is truncated by L. Below T*(L), the vari- 
ation of no comes mainly from the temperature depen- 
dence of the anomalous exponent, which is not as strong 
as the variation due to the exponential behavior of the 
correlation length. Consequently, dno j AT should ex- 
hibit a sharp minimum at the size-dependent temper- 
ature T* (L). Moreover, in a finite system, no cannot 
grow indefinitely as the temperature is lowered. With 
decreasing temperature (T — > 0), no must approach its 
(finite) T = value, which implies that dno/dT — > 0. 

Figure[3ta) depicts the derivative of the no for different 
system sizes vs T. The divergence of dno/dT is apparent. 
A sharp minimum develops and its location T* (L) ap- 
proaches T c as the system size increases. This is expected 
from the finite-size relation in Eq. (J7J . The scaling of the 
height of this minimum is studied in Fig. [3jb) , where we 
plot the absolute value of dno / dT\T* (l) vs L. The data 
follows the scaling relation in Eq. jSj), as made evident 
by a fit to the function g(L) = ao + aiL 7 / 4 In 3 (a2-L). In 
the inset in Fig. [3th), we show the finite-size scaling of 
T* (L). We observe that T* (L) is consistent with the 
scaling relation in Eq. (|7)l. which we use to obtain the 
critical temperature in the thermodynamic limit. We find 




FIG. 3. (Color online) (a) Derivative of the zero-momentum 
occupation no with respect to the temperature for different 
values of L. The inset shows no vs T. (b) Finite-size scaling 
of the height of the negative peak in dno/dT. The continuous 
line is a fit to the function g(L) = ao + aiL 7//4 In 3 (a2-I/). The 
inset shows the finite-size scaling of T* (L). 



T c /t = 0.701 ± 0.007. This value is compatible with the 
one found by performing the finite-size scaling of the su- 
perfluid stiffness. While this approach is obviously less 
accurate than the one discussed before for p Sl among 
other things because a numerical derivative is involved, 
the fact that it works extremely well is very important 
for current trapped ultracold gas experiments where the 
superfluid density cannot be measured. 

We note at this point that in the determination of 
Eq. ((5J), we have neglected multiplicative logarithmic cor- 
rections that affect the behavior of the zero-momentum 
occupation and thus its derivative with respect to the 
temperature (6?J, [jFj]. In fact, the exponent of the loga- 
rithm in Eq. ([5]) gets modified to 



dno 
~dT 



T*{L) 



L 7 / 4 m 3 ~ 2r L 
b 2 



(9) 



with r = -1/16 (Ref.lzFJ). However, this correction does 
not affect the determination of the critical temperature, 
which is based on the location of the position of the peak 
in the numerical derivative and the scaling relation in 
Eq. (J7J. Furthermore, the correction to the exponent of 
the logarithm is very small and, at least within the pre- 
cision of our simulations, its effect is hardly detectable. 



B. Three dimensions 

In order to determine the 3D phase diagram, we follow 
the same procedure as in 2D. In 3D, however, the super- 
fluid to normal transition belongs to the 3D XY univer- 
sality class. This transition, for the model in Eq. ([T]), has 
also been studied using QMC simulations in the past. T c 
for BEC was evaluated as a function of the density in 
Ref. [7l|. The onset of magnetization as a function of the 
magnetic field (or, in the bosonic language, the density 
as a function of the chemical potential) was investigated 
in Ref. [72|. Furthermore, the fate of the superfluid phase 
under the effect of an additional ring-exchange term was 
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FIG. 4. (Color online) (a) Superfluid stiffness of the 3D sys- 
tem for /j, — and several values of L. (b) Data collapse 
according to the relation in Eq. (|12l) . The inset shows the 
rescaled superfluid stiffness as a function of T. 



0.4 
0.3 
h§ 0.2 
0.1 








(a) 




L=8 






i L=12 


Q 




KJ L=16 






% L=20 






^\ L=32 






lH L=64 











1.6 



T/t 



1.5 




10 20 30 40 



FIG. 5. (Color online) (a) Condensate fraction in 3D for p, = 
and several values of L. (b) Data collapse according to the 
relation in Eq. (|14[) . The inset shows the rescaled condensate 
fraction as a function of T. 



studied in Ref. [73- Here, we determine the full phase 
diagram (shown in Fig. [1} as a function of the tempera- 
ture and the chemical potential. We begin by considering 
measurements of the superfluid stiffness. In d > 2 di- 
mensions, as the critical temperature is approached, the 
superfluid stiffness vanishes continuously as [z3| 

Ps ~ \T c -T\^ d - 2 >, (10) 

where the exponent v determines how the correlation 
length diverges when approaching the critical temper- 
ature, i.e., 

£~\T-T c \-». (11) 

As a result, at the critical temperature, the superfluid 
stiffness scales with the linear size of the system as 
p s L 2 ~ d . This, in turn, allows one to write the scaling 
hypothesis for the superfluid stiffness as a function of the 
system size and the temperature as 

p s L d - 2 = f(\T-T c \L 1 '^ , (12) 

which we utilize to determine the critical temperature. 
In Fig. 0|a) , we show results for the superfluid stiffness 
in a 3D lattice vs T for different system sizes. 

We numerically extract the scaling function F by 
studying the rescaled superfluid stiffness [left-hand side 
in Eq. ([T2])] vs the rescaled temperature (T - T C )L 1 / U . 
Classical Monte Carlo simulations yield the correlation 
length exponent v = 0.6717 ± 0.0001 [H, and v = 
0.6717± 0.0003 (HI, which we use to produce the collapse 
presented in Fig. QJb). With v at hand, it is enough to 
fix T c such that the best collapse of the data is achieved. 
Furthermore, the inset shows the rescaled superfluid stiff- 
ness as a function of temperature, which becomes system- 
size independent at the critical temperature, as implied 
by the scaling hypothesis in Eq. (fT2")) . Our best esti- 
mation of the critical temperature for p = is T c /t = 
2.0169 ± 0.0005 (to be compared with T c /t = 1.94 from 
Ref. [n| and more recently with T c /t = 2.016 ± 0.004 
from Ref. I75T ). We perform a similar analysis for dif- 
ferent values of the chemical potential to complete the 
three-dimensional phase diagram in Fig. [TJ 



Additionally, since the superfluid to normal phase 
transition in our model in 3D is accompanied by the 
emergence of true long-range order, one can study the 
transition by computing the condensate fraction fo as- 
sociated with the appearance of BEC. Following Pen- 
rose and Onsager [76J , the condensate fraction is defined 
as the ratio of the largest eigenvalue of the one-body 
density matrix to the total number of particles Nf,. For 
the system under consideration, condensation occurs to 
the zero-momentum state due to translational invariance, 
thus the condensate fraction is fo = uq/Ni,. The behav- 
ior of no can be obtained from the Fourier transform of 
the one-body density matrix in the long-distance limit, 
which in 3D is given by 

(a\a i+r )^r-^e W (-r/0. (13) 

Here, rj is the correlation function exponent, also known 
as the anomalous scaling dimension. On approach to T c , 
no diverges with the correlation length as [29[ 

no ~ e~ n - (14) 

In a finite system, this relation implies that the con- 
densate fraction vanishes at the critical point as fo ~ 
£-(!+>?) ; w hich we adopt to formulate the following scal- 
ing hypothesis for the condensate fraction 

foL 1+ ^ = F(\T-T c \L 1 l»y (15) 

In the determination of T c through the scaling relation 
in Eq. (IS]), we use the value n = 0.0381 ± 0.0002 [H2|. 
The results are summarized in Fig. [51 where a plot of 
the condensate fraction vs T is shown in panel (a). In 
Fig. [5jb), the data collapse of the rescaled condensate 
fraction foL 1+v vs the rescaled temperature is apparent. 
Furthermore, in the inset, one can observe that curves of 
the rescaled condensate fraction vs T become system-size 
independent at T c , as implied in Eq. (|15[) . This procedure 
results in a T c /t = 2.0167 ± 0.0005 for p, = 0, which is in 
remarkably good agreement with our previous estimate 
using the superfluid stiffness. 
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FIG. 6. (Color online) (a) dno/dT for different values of L. 
The inset shows no vs T. (b) Finite-size scaling of the height 
of the negative peak in dno/dT. The continuous line is a fit 
to the function g(\nL) = ao + ailnL. The inset shows the 
finite-size scaling of T* (L). 



1. Critical value from dno/dT 

Similarly to the 2D case, driQ / dT diverges in the vicin- 
ity of the superfluid to normal phase transition. It di- 
verges with the correlation length as 



dn 
~dT 



(16) 



Also, as in 2D, in a finite 3D system at a temperature 
T* (L) close to T c , the role of the correlation length is 
taken over by L when £ > L. The characteristic temper- 
ature T* (L) is given by 



T* (L) = T C + c'/L 

where d is a non-universal factor, 
scales with the system size as 



l/u 



(17) 



At T* (L), dn /dT 



driQ 
~dT 



_J^2-T]+1/V 



(18) 



T*{L) 



Furthermore, in a finite system, dno/dT reaches its min- 
imum value at T = T* (L) because the divergence of the 
correlation length can no longer be sustained. This is 
expected from the behavior of no vs T, shown in the in- 
set in Fig. Eta), where uq is first seen to increase as the 
temperature is lowered and then to saturate as T —> 0. 
The changes observed dno/dT in that low temperature 
regime originate in the smooth dependence of the corre- 
lation function exponent on the temperature, as opposed 
to the fast change produced by the strong divergence of 
the correlation length. Hence, once again, dno/dT ex- 
hibits a sharp minimum at the size-dependent tempera- 
ture T* (L) in Eq. (IT7|) and then goes to zero. 

In Fig. E{a), we display results for dno/dT vs T for 
different system sizes. The divergence in the derivative, 
anticipated by Eqs. ([To]) and ([T8"]). is confirmed by the 
presence of sharp minima that grow with system size. 
The finite-size scaling of the height of the sharp minimum 
in Eq. (|T8| is presented in Fig. E{b), where we plot the 
logarithm of the maximum height of \dno/dT\ vs InL. 



According to Eq. (|T8]). such a plot should turn into a 
straight line with a slope given by m = 2 — r\ + \jv. A fit 
of our data to the function o(lnL) = ao + a\ lni, yields 
ai = 3.47±0.01. The scaling relation given by Eq. (fT51) is 
thus confirmed as our value of a\ is compatible with the 
exponents from Ref. l52l which yield m = 3.450. The size 
dependence of the position of the peaks anticipated in 
Eq. (|T7|) is verified in the inset of Fig. [B^b) . Within this 
procedure, we find that the critical temperature in the 
thermodynamic limit is T c /t = 2.012 ± 0.002, which is in 
relatively good agreement with the one obtained through 
the finite-size scaling of both the superfluid stiffness and 
the condensate fraction. 

We conclude this section by mentioning that in deter- 
mining the critical temperature, we have used the lead- 
ing scaling forms and sublcading corrections to scaling 
have been neglected. For the 3D XY universality class, 
such corrections have been reviewed in Ref. WA We note 
that in our calculations, there is an excellent collapse of 
the data, which suggests that the effects of the sublead- 
ing corrections to scaling are small. Furthermore, the 
most accurate results obtained for T c follow from com- 
pletely independent measurements, i.e., the superfluid 
and condensate fractions. They agree within the error 
bars, which further supports the relevance of the scaling 
relations used. 



C. Mean field 

To gain an understanding of the effects of quantum 
fluctuations in our systems, we have also calculated the 
mean-field phase diagram for this model. We utilize the 
standard decoupling of the kinetic energy term in the 
Hamiltonian in Eq. (J) [ll] 



d\dj ~ a\^j 



(19) 



where = (ttj) is the condensate order parameter, to 
be determined self-consistently. The angle brackets de- 
note the usual thermal average. The above mean-field 
decoupling allows one to write a mean-field Hamiltonian 
for Eq. flTJ as 



(20) 



(i,f) 



For homogeneous systems, i.e., Vo = 0, Eq. (f20|) can 
be recast in the following manner, 



'MF 



= -2dt$ (a 



;,t 



fin, 



(21) 



where /imf is the mean-field Hamiltonian per lattice site. 
Note that in this case the superfluid order parameter can 
be taken to be real. The corresponding partition function 
at finite inverse temperature /3 is 



2e 



cosh/3 



(2dt$y 



(22) 
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A self-consistency condition for the superfluid order 
parameter can be derived by noting that 



A 7 

— = Af3dt{a)Z. 



(23) 



Using the relation (f2"B"|). we arrive at the equation that 
determines the order parameter $, 



(2dt<f>y 



dt tanh / 



(2eft$) 



(24) 



which is valid whenever $ > 0. We solve Eq. (|2~4"|) numer- 
ically and determine the superfluid region, $ > 0, as a 
function of the temperature and the chemical potential. 
The phase boundaries are determined as the values of p, 
and T for which <!> -> 0. For p, = 0, Eq. (JM]) reduces to 



2$ = tanh/3 2di$, 



(25) 



which is the equation that determines the mean-field 
magnetization of the Ising model in the absence of a 
magnetic field. The critical temperature is, of course 
T c /td = 1, which is quite different from the results of 
our quantum Monte Carlo simulations in two and three 
dimensions. 



IV. TRAPPED SYSTEMS 

In experiments involving ultracold atoms, an addi- 
tional trapping potential is necessary to contain the gas. 
While a qualitative (and sometimes a reasonably good 
quantitative) description of the trapped system can be 
obtained within the local density approximation (LDA) 
from the properties of the homogeneous system, this ap- 
proximation may breakdown in regimes of interest. In 
particular, the latter occurs at criticality, where the cor- 
relation length diverges and deviations from the LDA 
description can be large (29j . Furthermore, as we explain 
below, in trapped 2D systems care needs to be taken 
with the application of the Mermin-Wagncr-Hohenbcrg 
theorem. Therefore, we focus our attention on those two 
aspects, namely, the possibility to have BEC the in the 
presence of an additional external confining potential in 
2D, and the study of criticality in 2D and 3D. 



A. Absence of BEC in interacting 2D systems 

We mentioned in the Sec. U that homogeneous 2D sys- 
tems are special because thermal fluctuations destroy any 
order at finite temperature. However, harmonically con- 
fined non-interacting bosons can undergo BEC at finite 
temperature [79j]. In this case, the arguments by Mer- 
min et at arc not violated because condensation does 
not occur to the zero-momentum state, but to a single- 
particle eigenstate of the trapped system. One can then 
wonder whether finite-temperature BEC persists in the 



presence of interactions. By following analogous argu- 
ments to those in Ref. [8(1 we show below that interac- 
tions do preclude the formation of a condensate in the 
Bose-Hubbard model in the presence of the trap. This 
is so because there is a close connection between the for- 
mation of a condensate and the macroscopic population 
of the zero-momentum occupation, which is forbidden in 
2D at finite temperature. 

Generally speaking, the emergence of BEC is estab- 
lished through the evaluation of the condensate fraction 
/o, which is defined as the ratio of the largest eigenvalue 
of the one-body density matrix n m to the total number 
of particles iV&, 



, _ n M 
/0 ~ ~Nb' 



(26) 



If after taking the appropriate thermodynamic limit /o 
remains finite, then the system exhibits BEC. Otherwise, 
if it becomes zero, there is no condensation [Zf|. 

Alternative forms of the criteria expressed through 
Eq. (j2"6")l can be useful when the system is not spatially 
uniform; they are based on the following inequality (76| : 



n 2 M < ^ n 2 



(27) 



where n a are the eigenvalues of the one-body density ma- 
trix pij. We define the quantity 



Pij\ 



(28) 



which is just a lattice version of its analogous quantity 
defined on the continuum in Rcf. 

IB 

It follows from 

Eqs. (HU and (J2SJ that 



f$ < M < fo- 



(29) 



Therefore, if A 2 remains finite in the thermodynamic 
limit, then the system exhibits BEC. A further criterion 
can be defined and it depends on the quantity 



A 1 =(N b L d ) 'YAp^I. 



hi 



(30) 



Notice that (AiN b / L d ^ 2 is the square of the mean value 



of the function |py | , while A 2 (N b / L d ) 2 is the mean value 
of |pij | 2 . Since the variance of the function \pij\ is cither 
positive or zero, it follows that 

A\ < A 2 . (31) 

Now, since is a positive-semidefinite Hermitian ma- 
trix, its elements satisfy (ZfJ [8l| 

\pij I < Jp—p- < - fa + Pjj ) < aN b /L d , (32) 

where aN b / L d is an upper bound of the local density pa . 
By summing over i and j in Eq. (|32[) and the square of 
it, we find a lower bound for A\, 



A 2 < a A x . 



(33) 
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As long as the local density pa remains finite through- 
out the whole system, a can be taken to be finite and 
independent of N\,jL d . This, in turn, implies that if 
A\ > 0, then BEC takes place; otherwise if A\ = 0, then 
no BEC occurs [76|. Notice that if p^ > 0, then A\ coin- 
cides with the ratio of the zero-momentum occupation to 
the total number of particles, i.e., the fraction of parti- 
cles in the system that condenses to the zero-momentum 
state. Since in two dimensions no/Nb vanishes because 
of the Mermin et al. theorem, then Ax is zero too. In the 
specific case of the Bose-Hubbard model in the presence 
of an inhomogeneous potential in thermal equilibrium, 
we have that > 0. Furthermore, the density is fi- 
nite everywhere across the system because of the on-site 
interaction, implying that A\ = 0. 

Hence, even in the presence of the trap, there is no 
condensation in the 2D Bose-Hubbard model at finite T . 
Note that this argument does not preclude condensation 
in the non-interacting limit, where the density can di- 
verge at the minimum of the inhomogeneous potential in 
the thermodynamic limit and BEC can indeed occur to 
the lowest single-particle eigenstate, but not to the zero- 
momentum state. Moreover, the criteria above implies 
that for the Bose-Hubbard model in d > 2 in thermal 
equilibrium, condensation to any state has to be accom- 
panied by condensation to the zero-momentum state. 

In our proof, we have stated that for the Bose-Hubbard 
model in thermal equilibrium p^ > holds. We now 
present two independent arguments for why pij > 0. The 
first one is based on the fact that the matrix elements of 
the von Neumann's statistical operator in the position 
representation are strictly positive (82j . Since the one- 
body density matrix corresponds to a partial trace of the 
von Neumann's statistical operator [76j |, it follows that 
its elements are positive too. A rather technical, but yet 
rigorous, argument is based on the series expansion rep- 
resentation of the one-body density matrix that we used 
in our Monte Carlo implementation. Within this rep- 
resentation, the measurements of the one-body density 
matrix are based on the extension of the configuration 
space where these off-diagonal quantities are well defined 
[59(. In such extended space, the one-body density ma- 
trix is represented as the sum of strictly positive matrix 
elements (hence pij > 0) which are, in turn, efficiently 
sampled during the construction of the loop operators in 
the directed- loop update algorithm (83j . 



B. Two dimensions 

1. Local compressibility 

Exactly as in the homogeneous system, even though 
there is no condensation in 2D, a supcrfluid phase is ex- 
pected in the trapped system at low temperatures. Be- 
cause of the inhomogencity introduced by the confining 
potential, the coexistence of space separated normal and 
supcrfluid domains can occur at intermediate tempera- 



ture. In that case, there must be a region in the trap 
where supcrfluidlike domains transition into normal ones. 
Within the LDA, this region is such that the local chem- 
ical potential pi coincides with the critical p of the bulk 
system for the normal-to-superfluid phase transition. 

Based on this idea, Zhou and collaborators proposed 
a method to identify the phase boundaries of the homo- 
geneous system from a high-resolution scan of the local 
density p(r) across the confined system (3l| . This method 
requires the determination of the local compressibility de- 
fined as 



Kdiff (r) 



1 dp (r) 
2V r dr : 



(34) 



and relies on the expectation that the local density profile 
p(r), as well as the local compressibility Kdiff (r), can be 
well approximated by their bulk values through the LDA. 
The existence of sharp features in the local compressibil- 
ity at specific locations in the trap is then associated with 
phase transitions occurring in the homogeneous system 
as a function of the chemical potential. This method is 
expected to l ie accurate in the limit of very shallow traps 
where the contribution from density gradients due to the 
trapping potential are small f&H - 

In Fig. [Jj we present QMC results for the density pro- 
file of a 2D trapped system, as well as the local com- 
pressibility, as a function of the distance from the center 
of the trap. The expected sharp features in the local 
compressibility due to critical fluctuations are smoothed 
by finite-size effects. They are replaced by a rounded 
maximum, which can be associated with the supcrfluid 
to normal transition [84j . The location of the maxi- 
mum r c is connected to the critical chemical potential 



through p c 



Vor^. For the case in Fig. [JJ we get 



p c /t = — 3.57± 0.03. This value is to be contrasted with 
p c /t = —3.5, which we obtained in the homogeneous sys- 
tem calculations. As T increases, however, the agree- 
ment between the estimates of the critical chemical po- 
tential based on the local compressibility and the results 
of the homogeneous system worsens. For instance, for 
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FIG. 7. (Color online) (a) Two-dimensional density at T/t = 
0.2012 and a trapping potential Vo/t — 0.0003, for p = in 
the center of the trap, (b) The corresponding density profile 
p(r), as well as the local compressibility Kdiff (r), as a function 
of the distance from the center of the trap r. All distances x, 
y, and r are measured in units of a while the local compress- 
ibility is measured in units of 1/t 
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FIG. 8. (Color online) (a) no as a function of T and /3 in a 
trapped 2D system with Vo/t — 0.00125, /i = in the center 
of the trap, and L — 128. (b) Derivatives of no with respect 
to P and with respect to T. 

T/t = 0.4562, wc find that fx c /t = 2.99±0.04, as opposed 
to the homogeneous system result where fx c /t = 2.5. This 
occurs presumably because, closer to the tip of the su- 
perfluid lobe, critical fluctuations are stronger, and thus 
larger violations of the LDA are expected. 

2. Momentum distribution function 

Another quantity that can be measured in experiments 
with ultracold atoms is the momentum distribution func- 
tion. At fixed chemical potential (/x < 0), when lowering 
T, the normal-to-supcrfluid crossover in the trapped sys- 
tem proceeds via the creation and growth of a supcrfluid 
domain in the center of the trap. (The rate of growth of 
the supcrfluid domain will depend on the functional form 
and strength of the confining potential.) Hence, the zero- 
momentum state becomes increasingly populated. As fol- 
lows from the discussion for finite homogeneous systems, 
it is expected that as T decreases and approaches T c for 
the normal-to-superfluid transition in the center of the 
trap, the rate of growth of no will increase. Below T c , 
on the other hand, dno/dT will eventually decrease be- 
cause of the finite extend of the system imposed by the 
confining potential. If T is lowered well below T c , then 
almost the entire system will become superfluid and the 
observablcs will saturate their (finite) zero-tempcraturc 
values. 

Hence, just as in the homogeneous case, one can at- 
tempt to estimate T c for the superfluid to normal phase 
transition for the density in the center of the trap by 
measuring the temperature at which the rate of change 
of no is extremal. This approach provides an accurate 
estimate for the homogeneous system and it is expected 
to be accurate in confined systems with shallow trap- 
ping potentials. Figure [H^a) depicts the evolution of tiq 
vs T as well as the inverse temperature /3 of a harmoni- 
cally confined 2D system with V /t = 0.0015 (L = 128) 
and fx = in the center of the trap. In Fig. EJb), we 
show dno/dT which, as expected, exhibits a minimum 
located at T/t = 0.66 ± 0.02. This temperature is com- 
patible with the value of T c /t obtained for the homoge- 



neous case where, after a finite-size scaling, we obtained 
T c /t = 0.685 ± 0.001. Our estimate derived from the 
study of a single trapped system is about 4% off the value 
of the homogeneous system. 

One can perform the same analysis based on measure- 
ments of no, but now as a function of the inverse tem- 
perature /?. In that case, one expects a maximum in the 
derivative dno/d/3 instead of a minimum. In general, for 
finite and not very large systems, the position of such 
maximum f3 c will not coincide with 1/T C obtained from 
the minimum of dno/dT. Overall, we find that for the 
system sizes available to our QMC simulations, the anal- 
ysis based on dno / d/3 provides more accurate estimates of 
the critical temperature than the one based on dno/dT. 
Furthermore, the maximum found in dno/d/3 is consis- 
tently sharper and better defined with respect to the 
minimum found for dno/dT which instead is shallower 
and broader, and thus harder to detect and numerically 
less reliable. 

Based on measurements of dno /dp presented in Fig. 
Etb) on the same system with V /t = 0.0015 (L = 128), 
u = 0, we find T c /t = 0.72 ±0.02, which is also very close 
to the critical temperature of the homogeneous system. 
When the maximum is sharply defined, in the limit of 
very shallow traps with large numbers of bosons, the two 
approaches are expected to coincide (i.e., their difference 
is due to finite-size effects). As a matter of fact, for the 
homogeneous 2D and 3D systems in Sec. IIII1 where the 
minima of dno/dT are sharp, we find that the analysis 
using dno/dT and dno/d/3 yields essentially the same re- 
sults for T c . In the Appendix, we provide an analytic un- 
derstanding of this in terms of a simple function. There- 
fore, for the determination of the phase diagram based 
on measurements in harmonically confined systems, we 
consider only measurements based on dno/d/3. 

In Fig. [9l we summarize our results for the determina- 
tion of the critical parameters with the local compressibil- 
ity as well as with the derivative of the zero-momentum 
occupation with respect to f3, and contrast them with 
the phase diagram of the homogeneous system. Clearly, 
all methods work well for large values of fi/t and small 
values of T c /t (equivalent to approaching the continuum 
limit in a lattice system). Close to the tip of the su- 
pcrfluid region, the method based on no performs much 
better than the one based on Kdiff (r). 

At the tip of the superfluid lobe, where the size ef- 
fects are expected to be the strongest, we observe that 
as the size of the system is increased (or the strength 
of the trap is decreased) , keeping constant the chemical 
potential in the center of the trap, the estimate of the 
critical temperature decreases approaching the result in 
homogeneous systems. 



C. Three dimensions 

We now turn our attention to the study of criticality 
in 3D trapped systems. We make use of the same ideas 
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FIG. 9. (Color online) Estimate of the critical points based 
on the local compressibility (black dots based on a system 
with L — 256) as well as the derivative of the zero-momentum 
occupation with respect to /3 (red diamonds based on a system 
with L = 128). At the tip of the superfluid lobe we include 
further results for different system sizes and trap strengths 
(yellow pentagon L = 32, blue triangle L — 64, violet empty 
circle L — 256). The phase diagram of the homogeneous 
system is also shown. 

developed for 2D system to extract the critical param- 
eters, i.e., measurements based on the zero- momentum 
occupation as well as on the local compressibility. 

Additionally, in 3D, we can utilize a method that is 
based on the analysis of the shape of the central peak for 
the momentum distribution. With it, one can construct 
aquantity that exhibits a minimum at the critical point 
[23 ]. The idea behind this method is that close to crit- 
icality, the momentum distribution develops a bimodal 
structure whose evolution as a function of temperature 
contains information about the formation of a superfluid 
region in the center of the trap. At T c , when a superfluid 
domain begins to form, the major contribution to the oc- 
cupation of the zero- momentum state comes from regions 
that are not critical, i.e., from regions that are far away 
from the center of the trap. However, the derivatives of 
the momentum distribution d m rik/dk m are critical, in the 
sense that they can be understood in terms of an LDA 
integral that diverges at the center of the trap, where 
the system is critical. Based on that idea, the follow- 
ing quantity was devised in order to extract the critical 
temperature (29j : 

Q (T) = (n - n fcm „) (k niax ) s , (35) 

where /c max is the momentum at which \drik/dk\ is max- 
imum and the exponent s > 2 — n. In Ref. |29| . it was 
shown that Q (T) should exhibit a minimum at the crit- 
ical temperature T c . 

We plot Q (T) vs T in Fig. ITOTal Q (T) exhibits a 
minimum at T c /t = 2.04 ± 0.03. In the inset, we show 
the evolution of the momentum distribution function as 
the temperature of the system is reduced. This result is 
compatible with the critical temperature found for the 
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FIG. 10. (Color online) (a) The quantity Q (T) as a function 
of temperature extracted from the momentum distributions 
shown in the inset. The exponent in Eq. 1)351) has been set 
to s = 3. (b) Derivatives of no with respect to f3 and with 
respect to T. The three-dimensional system is prepared with 
Vo/t = 0.04, fi = 0, and L = 32. 



homogeneous system T c /t = 2.0169 ± 0.0005. In princi- 
ple, similar ideas as the ones presented in Ref. [29J could 
be used to devise a quantity Q (T) to locate the critical 
parameters in 2D. In that case, however, the structure of 
the momentum distribution is different because the tran- 
sition is in another universality class. As a result, the 
LDA integrals for the central peak and the derivatives of 
the momentum distribution get substantially modified. 
We find that both the central peak and the derivatives of 
rife are critical in 2D because the LDA integrals of those 
quantities diverge in the center of the trap where the sys- 
tem is critical. Hence, one cannot define a Q (T), as done 
in 3D, that will exhibit a minimum at T c . 

In Fig. nUT b). we also display results obtained for 
dno/d(3 (dno/dT) in the same system. The tempera- 
ture at which the maximum (minimum) occurs for those 
quantities exhibits a larger deviation from T c , from the 
homogeneous case, than Q(T). However, with increasing 
system size, we find that the maxima of dno/d(3 (minima 
of dno/dT) slowly approach the homogeneous result. In 
experiments where the system sizes are much larger than 
the ones studies here, we expect that dno/dT and dno/d(3 
will both produce accurate results for T c . 

In Fig. 11 11 we present a summary of our estimates of 
the critical parameters based on the local compressibil- 
ity, the derivatives of no with respect to (3, and Eq. (|3"5j) . 
The method based on Q (T) is found to be more accurate 
than those based on dno /dj5 and the local compressibility. 
This is understandable because the former approach uses 
precise information of the nature and universality class 
of the transition in 3D. Nevertheless, as argued before, 
we anticipate that if one decreases the strength of the 
confining potential and increases the number of bosons, 
as to reach the system sizes that are studied experimen- 
tally, then dno/d(3 will provide accurate results (at least 
similar to the ones obtained in 2D). This effect is stud- 
ied in Fig. 1111 where we show the evolution of the critical 
temperature at the tip of the lobe as a function of sys- 
tem size. As the strength of the confining potential is 
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FIG. 11. (Color online) Estimates of the critical points based 
on the local compressibility (black dots; based on a system 
with L = 64) as well as the derivatives of the zero-momentum 
occupation with respect to j3 (red diamonds based on a system 
with L = 32). At the tip of the superfluid lobe, we include 
further results for different system sizes and trap strengths 
(yellow pentagon: L = 16; blue triangle: L = 20). The esti- 
mates based on Eq. (|35p are shown by purple empty triangles. 
The phase diagram of the homogeneous system is also drawn 
with green empty squares. 

decreased and the size of the system is increased, the es- 
timate of the critical temperature based on dno / df3 tends 
to increase and approach T c in the homogeneous system. 
The method based on the local compressibility is found 
to be inadequate close to the tip of the lobe. This is be- 
cause the maximum of Kdiff (?") becomes very broad and 
finite-size effects are stronger. In that regime, one also 
needs a higher accuracy in the determination of the den- 
sity in order to accurately compute the local compress- 
ibility In spite of this, in 3D, the method based on the 
local compressibility yields more accurate results than in 
2D (compare Figs. l9l and [TTj) . 



more, we computed the phase diagram using mean-field 
theory and found it to be quantitatively quite different 
from the results of numerically exact QMC simulations 
in 2D and 3D. Hence, for this model, thermal and quan- 
tum fluctuations are strong even in three dimensions, and 
mean-field theory is a poor approximation. 

In the presence of an additional confining potential, 
we proved that the Bose-Hubbard model does not exhibit 
finite-temperature BEC in two dimensions, provided that 
density remains finite across the entire system in the 
thermodynamic limit. Moreover, we considered measure- 
ments of the critical temperature and chemical poten- 
tial of the homogeneous system based on experimentally 
measurable quantities such as the momentum distribu- 
tion function and the local density profile. The accuracy 
of each method discussed depends on the dimensionality 
of the system and the range of temperatures and chem- 
ical potentials considered. In two dimensions, we found 
that the approach introduced in this work, based on the 
derivatives of no with respect to /3, is accurate in all re- 
gions of the phase diagram. A method based on the 
measurement of the local density was found to be reliable 
when T c is low, while close to the tip of the superfluid 
lobe this approach is less effective even when the trap is 
very shallow. This can be understood to be due to the 
strong deviations from the LDA close to the tip of the 
superfluid lobe. A quantitative account of these devia- 
tions based on trapped finite-size scaling, as presented 
in Ref . [85| and [8fj| . would in principle allow one to per- 
form an accurate size-scaling analysis in the presence of 
the confining potential, which might potentially improve 
the capabilities of the methods based on the measure- 
ments of the density profile. The accuracy of the latter 
method improves in 3D, but still remains inadequate as 
one approaches the tip of the superfluid lobe. In three 
dimensions, the approach based on Q (T) was found to 
be the most accurate. 



V. CONCLUSIONS 

We have presented a detailed study of the finite tem- 
perature phase diagram of strongly correlated bosons in 
the hard-core limit (or the XY model) in two and three 
dimensions. The critical parameters in the homogeneous 
case were determined through a finite-size scaling analy- 
sis of the superfluid stiffness and the condensate fraction. 
We introduced an approach to estimate the critical tem- 
perature from measurements of no in finite systems. It 
makes use of the behavior of the derivative dno / dT and 
we derived finite-size scaling relations that can be used to 
extrapolate the results to the thermodynamic limit. This 
approach can be applied to systems that exhibit a diverg- 
ing zero-momentum occupation in any dimension, irre- 
spective of the universality class to which the transition 
belongs. We showed that this method is also accurate 
in 2D, where the system does not exhibit BEC. Furthcr- 
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Appendix A: Differences between dno/dT and dno/d/3 

We briefly illustrate, by means of a simple analysis, 
why the estimate of the critical temperature based on 
dno/dT differs from the estimate based on dno/d(3. We 
also discuss under which conditions the two estimates 
should approach each other. 

We consider no(/3) to be the zero momentum occu- 
pation in the vicinity of /3 C . Its first derivative, which 
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exhibits a maximum at j3* , can be written as 



1//3* given by 



driQ 
~dj3 



(Al) 



where the curvature of the parabola is a < 0, the height of 
the maximum is do , and /3 is assumed to be very close to 
(3*. If instead we now compute dn$/dT, we anticipate a 
minimum of this function located at a temperature T* =/= 



1 



3a/3* 



8d Q 



■ia 



(A2) 



In general, the position of the minimum as a function 
of T depends on the position of the maximum /?*, its 
curvature a, and its height do. However, in the limit 
of very large system sizes and very shallow traps, one 
expects the maximum of the derivative dn /d/3 to be very 
sharp. In our simple example, this regime corresponds to 
a large value of the curvature, i.e., |do/ a l "C /3* 2 , which 
implies that T* ~ 
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